arXivil504.06619vl [astro-ph.CO] 24 Apr 2015 


Mon. Not. R. Astron. Soc. 000, 000-000 (0000) Printed April 28, 2015 (MN style file v2.2) 


The BlueTides Simulation: First Galaxies and Reionization 

Yu Feng*^’^, Tiziana Di-Matteo^, Rupert A. Croft^, Simeon Bird^, Nicholas Battaglia^’^ 
Stephen Wilkins^ 

^McWilliams Center for Cosmology, Carnegie Mellon University, Pittsburgh PA, 15213 
“^Berkeley Center for Cosmologieal Physies, University of California, Berkeley, Berkeley CA, 94720 
^Department of Astrophysieal Seienees, Prineeton University, Prineeton, NJ 08544 

Astronomy Center, Department of Physies and Astronomy, University of Sussex, Brighton, BN19QH, UK 


April 28, 2015 


ABSTRACT 

We introduce the BlueTides simulation and report initial results for the luminosity 
functions of the first galaxies and AGN, and their contribution to reionization. Blue¬ 
Tides was run on the Blue Waters cluster at NCSA from z = 99 to z = S.O and 
includes 2x7040^ particles in a 400 h~^Mpc per side box, making it the largest hydro- 
dynamic simulation ever performed at high redshift. BlueTides includes a pressure- 
entropy formulation of smoothed particle hydrodynamics, gas cooling, star formation 
(including molecular hydrogen), black hole growth and models for stellar and AGN 
feedback processes. The star formation rate density in the simulation is a good match 
to current observational data at z ^ 8 — 10. We find good agreement between observa¬ 
tions and the predicted galaxy luminosity function in the currently observable range 
—18 < Mjjv < —22.5 with some dust extinction required to match the abundance 
of brighter objects. BlueTides implements a patchy reionization model that produces 
a fluctuating UV background. BlueTides predicts number counts for galaxies fainter 
than current observational limits which are consistent with extrapolating the faint 
end slope of the luminosity function with a power law index 1.8at^^8 

and redshift dependence of a ^ (1 -h The AGN population has a luminosity 

function well fit by a power law with a slope a ^ —2.4 that compares favourably 
with the deepest CANDELS-Goods fields. We investigate how these luminosity func¬ 
tions affect the progress of reionization, and find that a high Lyman-a escape fraction 
(/esc ^0.5) is required if galaxies dominate the ionising photon budget during reion¬ 
ization. Smaller galaxy escape fractions imply a large contribution from faint AGN 
(down to Muv = —12) which results in a rapid reionization, disfavoured by current 
observations. 
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1 INTRODUCTION 


Recent deep observations using the Hubble Space Telescope 
have detected a plethora of objects at ever higher redshift 
(Bouwens et al. 2014a| 


Oesch et al. 2014a McLeod et al. 


2014) and measured the galaxy UV luminosity function at 
z < 10. At these redshifts microwave background measure¬ 
ments suggest that a substantial fraction of the Universe 


is still neutral (Hinshaw et al. 2013), and these observa¬ 


tions may therefore probe the epoch of reionization. In the 
near future, next-generation space missions such as JWST 
( Gardner et al.|20Q6 ) and WFIRST ( Spergel et al.|2013 ) will 
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increase the number of available samples by several orders of 
magnitude. These are expected to detect the sources which 
produce the ionizing photons that drive reionization. 


The formation of these objects is driven by non-linear 
gravitational collapse and so understanding them requires 
cosmological hydrodynamic simulations. The presence of an 


ionizing background may affect galaxy formation (Madau, 


Pozzetti & Dickinson 1998). To model the processes gov¬ 


erning reionization it is desirable to simultaneously include 
scales of a few hundred Mpc, the characteristic size of ioniza¬ 
tion bubbles, and, to resolve the formation of galactic halos, 
scales of a few kpc. 


We present results from BlueTides, the largest cosmo¬ 
logical hydrodynamic simulation yet performed, enclosing 
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a box 400/i“^Mpc on a side, with a smoothing length of 
1.5 h~^Kpc, and including 2 x 7040^ particles- a total of 0.7 
trillion particles. Our high resolution allows us to study the 
formation of disc galaxies ( Feng et al.|2015a ), while the large 
volume allows study of the progress of reionization. 

We include a number of physically relevant processes, 
including star formation incorporating the effects of molec¬ 


ular hydrogen formation (Krumholz & Gnedin 2011), en¬ 


ergetic feedback from supernovae (Okamoto et al. 2010), 
and feedback from super-massive black holes (Di Matteo, 
Springel Hernquist|2005 ). We include for the first time in a 
simulation of this size a model for patchy reionization which 


varies the optical depth based on local density (Battaglia 
et al.|2013| |. 


Several large volume simulations have recently been 
performed. In particular, MassiveBlack I was the largest vol¬ 
ume hydrodynamic simulation Di Matteo et al. (2012) to 


study reionization at z > 4.75. MassiveBlack II simulated 
a 100/i“^Mpc volume at substantially improved resolution 
to z = 0 ( Khandai et al.|[2014 ). Illustris included a volume 
and resolution comparable to MassiveBlack II, but with im¬ 
proved prescriptions for the effect of energy injection from 
supernovae ( Okamoto et al.|[2010 ). The Eagle simulation is 
similar in size to MassiveBlack II and Illustris, but with a dif¬ 
ferent approach to sub-grid modeling which allows improved 
agreement with observations by weakening requirements for 
numerical convergence with resolution ( Schaye et al.||2015 ). 
Concurrently, dark matter only simulations have continued 
to increase in size, reaching loads of trillions of particles 


Habib et al. (2014). 


BlueTides is based on the simulation code used in Mas¬ 
siveBlack I & II, P-Gadget3 ( Springel|2005 Di Matteo et al. 
2012 Khandai et al. 2014). The simulation encloses a volume 


comparable to MassiveBlack I, has a resolution comparable 
to MassiveBlack II, and includes a stellar feedback model 
similar to that of Illustris. Our particle load is ten times 
larger than that of MassiveBlack I, which was previously 
the largest hydrodynamic simulation. 

In Sectionwe present our methods, explaining briefly 
the computational techniques necessary to perform a sim¬ 
ulation of this magnitude. In Sections and ^ we examine 
the basic statistics of objects within the simulation. We show 
the galaxy and AGN UV luminosity functions and star for¬ 
mation rates, and we present fits to these functions for easy 
comparison to future observations. In Sectionwe compute 
the sources of ionising photons within our model and use 
them to examine features of reionization. Finally we con¬ 
clude in Section |6] 


2 BLUETIDES: SOFTWARE AND SUB-GRID 
PHYSICS 

The BlueTides simulation was performed on the Blue Waters 
cluster at the National Genter for Super-computing Appli¬ 
cations (NGSA). We operated the production run on a total 
of 648,000 Gray XE compute core of Blue Waters. This is the 
largest cosmological hydrodynamic simulation to date, con¬ 
taining a simulation volume roughly 300 times larger than 


the largest observational survey at redshift 8—10 (Trent! 
et al.|2011 |. This extraordinary size allows us to easily com¬ 
pare our results to current and future observations, and ob¬ 


tain a representative sample of the first galaxies which may 
have driven reionisation. A visual overview of the simulation 
is shown in Figure 

Halos in BlueTides are identified using a Friends-of- 
Friends algorithm with a linking length of 0.2 times the mean 


particle separation (Davis et al. 1985). We have not per¬ 


formed sub-halo or spherical over-density finder algorithms 
on BlueTides due to limits on the scalability of current im¬ 
plementations. We will however investigate the mass func¬ 
tion with spherical over-density finders in a follow up work. 


2.1 Computing: Improved performance at 
Peta-Scale 

At the particle numbers reached by our simulation, any un¬ 
necessary communication overhead can easily become a sig¬ 
nificant scalability bottleneck. In order to allow the simula¬ 
tion to fully utilize the available computational power, we 
implemented several improvements to the speed and scala¬ 
bility of the code. Here we briefly list the substance of the 
most important changes, deferring a fuller description to 
Fengl (|20I5|. 


First, we substantially improved the scalability of the 
threaded tree implementation, which computes short-range 
particle interactions. This includes the gravitational force on 
small scales, hydrodynamic force, and the various feedback 
processes. Our improved routine scales to > 30 threads and 
at 8 threads is twice as fast as the previous implementa¬ 
tion in P-Gadget3. Scalability improvements were achieved 
by eliminating OpenMP critical sections in favour of per- 
particle POSIX spin-locks, essentially making the thread ex¬ 
ecution wait-free at high probability. 

Second, we replaced the default particle mesh grav¬ 
ity solver (used to compute long-range gravitational forces) 
based on FFTW with a gravity solver based on a 2-d tile 
FFT library, PFFT ( Pippig|[20I3 ). This allows a more effi¬ 
cient decomposition of particles to different processors, sig¬ 
nificantly improving the load and communication balance. 
To further simplify the book-keeping and reduce memory 
overhead, we switched to using Fourier-space finite differenc¬ 
ing of gravity forces, as used by the N-body gravity solver 
HACC dHabib et al.1|M4 |. We also reduced the communi¬ 
cation load by applying a sparse matrix compression of the 
local particle mesh. In concert these changes removed the 
PM step as a scalability bottleneck. 

Third, at the high redshift covered by BlueTides, only 
a small fraction of the particles are in collapsed halos. Thus, 
to ease analysis, we produce two digest data sets on the fly 
in addition to the full simulation snapshots: I) Particles- 
In-Group (PIG) files, which contain the attributes of all 
particles in the over-dense regions detected by the Friend 
of Friend groups; 2) Subsample files, containing a fair sub¬ 
sample of I/I024 of the dark matter and gas particles, and 
a full set of star and black hole particles. 

Fourth, we implemented a histogram based sorting rou¬ 
tine to replace the merge sort routine originally present in P- 
Gadget3. Histogram based sorting routines have been shown 
to perform substantially better at scale (Solomonik & Kale 


2010), a result confirmed by our experience in BlueTides 


(Feng et al. 2015b). Particles must be sorted when construct¬ 


ing the FoF group catalogues, and our new sort routine sped 
up FoF catalogue generation times by a factor of 10. The 
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Figure 1. Overview of the BlueTides Simulation. Left inset: a field of view from the BlueTides simulation at 2 : = 8 with the same total 
area as the BoRG survey ( [Trent i et al.||2011| |. Right inset: the field of view of the entire BlueTides simulation at 2 ; = 8. Background: 
galaxies in BlueTides. Even rows: top, front, and side views of the gas component of FoF halos. Odd rows: top, front, and side views of 
the star formation rate surface density of FoF halos. 


source code of this sorting routine, MP-Sort, is available from 
http://github.com/rainwoodinan/MP-sort to facilitate in¬ 
dependently reuse. 

Finally, we implemented a new snapshot format (Big- 
File) that supports transparent file-level striping and sub¬ 
stantially eases post-production data analysis compared to 
multiple plain HDF5 files. Until recently, file-level striping 
(bypassing MPIO) has been the only way to achieve the 
full 10 capability of Lustre file system for problems at our 
scale. We release the library for accessing the BlueTides sim¬ 
ulation at http://github.com/rainwoodman/bigfile, to¬ 
gether with Python language bindings for post-simulation 
data analysis. 


2.2 Physics: Hydrodynamics and sub-grid 
modelling 

BlueTides uses the hydrodynamics implementation de¬ 


scribed in Feng et al. (2014). We adopt the pressure-entropy 


formulation of smoothed particle hydrodynamics (pSPH) to 
solve the Euler equations ( Hopkin^|2013 Read, Hay field 
AgertzpOlO ). The density estimator uses a quintic density 
kernel to reduce noise in SPH density and gradient estima¬ 


tion (Price 2012) 


i^pOl 

4e (3 li 


Table [U lists the basic parameters of the simulation. 


Initial conditions are generated at 2 : = 99 using an initial 
power spectrum from CAMB (Lewis & Bridle 2002). Star 


formation is implemented based on the multi-phase star for¬ 


mation model in Springel & Hernquist (2003), and incor¬ 
porating several effects following Vogelsberger et ar|(2013|. 
Gas is allowed to cool both radiatively following jKatz, Wein-j 
berg & Hernquist (1996) and via metal cooling. We approx¬ 
imate the metal cooling rate by scaling a solar metallicity 
template according to the metallicity of gas particles, fol¬ 


lowing Vogelsberger et al. (2014). We model the formation 


of molecular hydrogen, and its effect on star formation at 
low metallicities, according to the prescription by jKrumholzj 
& Gnedin (2011). We self-consistently estimate the fraction 
of molecular hydrogen gas from the baryon column density, 
which in turn couples the density gradient into the star- 
formation rate. 

A SNII wind feedback model ( Okamoto et al.pOlO ) is 
included, which assumes wind speeds proportional to the 
local one dimensional dark matter velocity dispersion ctdm: 


Vw — TiwCTDM , 


( 1 ) 


where Vw is the wind speed. Kw is a dimensionless parameter, 
which we take to be 3.7 following Vogelsberger et aL] ( |2013| ). 

We model feedback from active galactic nuclei (AGN) 
in the same way as in the MassiveBlack I V II simulations. 
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Name 

Value 

Notes 

Ha 

0.7186 


HlVatter 

0.2814 

Baryons + Dark Matter 

H Bar yon 

0.0464 


h 

0.697 

Hubble parameter in units of lOOkm/s/Mpc 

crs 

0.820 


ris 

0.971 


Ebox 

400 /i“^Mpc 

Length of one side of the simulation box. 

Nparticle 

2 X 7040^ 

Total number of gas and dark matter particles in the initial conditions 

Mdm 

1.2 X lO'^/i-^M© 

Mass of a dark matter particle. 

Mgas 

2.36 X lO^/i-^M© 

Mass of a gas particle in the initial conditions. 

m 

I.O 

SPH smoothing length in units of the local particle separation^ 

^grav 

1.5 h~^Kpc 

Gravitational softening length. 

-^Generation 

4 

Mass of a star particle as a fraction of the initial mass of a gas particle. 

egy^/egyo 

I.O 

Fraction of supernova energy deposited as feedback. 

K,w 

3.7 

Wind speed as a factor of the local dark matter velocity dispersion. 

FsnII,51 

I.O 

Supernova energy in units of 10^^ erg/s 

<u 

5 X 10^ /i-^M© 

Seed mass of black holes. 

Kalo 

5 X 10^2/i-iM© 

Minimum halo mass considered in black hole seeding. 

VBH 

0.05 

Black hole feedback efficiency. 


t A value of 1.0 translates to 113 neighbour particles with the quintic kernel used in BlueTides. 


Table 1. Parameters of the BlueTides Simulation 


using the super-massive black hole model developed in 


Matteo et ah (2005). Super-massive black holes are seeded 
with an initial mass of 5 x 10^ h~^MQ in halos more mas¬ 
sive than 5 x 10^° while their feedback energy is de¬ 

posited in a sphere of twice the radius of the SPH smoothing 
kernel of the black hole. 

The large volume of BlueTides allowed us to include 
some of the effects of “patchy “ reionization, where the am¬ 
plitude of the ultraviolet background is spatially variable. 
We model patchy reionization using a semi-analytic method 
based on hydrodynamic simulations performed with radia¬ 
tive transfer (for more details see Battaglia et al.|2013 ). This 
method uses an evolved density field calculated from the 
initial conditions using second order Lagrangian perturba¬ 
tion theory to predict the redshift at which a given spatial 
region will reionize. In our fiducial reionization model, we 
set the mean reionization redshift at z ~ 10 based on the 
measurement of the optical depth, r, from the WMAP 9 
year data release ( Hinshaw et al.|2013 ). In regions that have 
been reionized, we assume the UV background estimated by 


Faucher-Giguere et al. (2009). The global neutral fraction 


in BlueTides evolves smoothly as a function of redshift, as 
seen in Figure 


3 STAR FORMATION RATE 

Figure shows the global star formation rate in BlueTides, 
together with observational constraints and several other 
simulations. We show the total star formation density in the 
whole volume (dotted lines) and the total star formation rate 
for halos with SFR > 0.7M©/year (solid line). The latter 
is directly comparable with current observations and corre¬ 
sponds to the current observational limit of Muv = —18. 
The SFR density in BlueTides smoothly increases with de¬ 
creasing redshift. 

Several Hubble ultra deep field surveys have given esti¬ 
mates on the star formation rate density due to halos with 



Figure 2. Global neutral hydrogen fraction as a function of red¬ 
shift in BlueTides. 


Muv < —18. BlueTides halos at the observational limit typ¬ 
ically contain a few thousand particles, and are thus well 
resolved. The BlueTides predictions for the star formation 
rate is in good agreement with current observations, with 
the caveat that at these high redshifts observational uncer¬ 
tainty remains high. 

Figure also compares the SFR density in BlueTides to 
that from two other recent (albeit smaller volume) simula¬ 
tions: Illustris ( |Genel et al.||2014 Vogelsberger et ^|2014 ), 
and MassiveBlack II (Khandai et al. 2014| . Illustris uses 
similar prescriptions for sub-grid feedback, but a different 
solver for the Euler equations, while MassiveBlack II uses 
substantially different sub-grid modelling, which is less effec¬ 
tive at suppressing star formation in faint objects. Neither of 
the other simulations include patchy reionization. The star 
formation rates for BlueTides and Illustris agree very well. 
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Figure 3. The Star Formation Rate Density in the BlueTides simulation. Solid blue: star formation rate density of halos with star 
formation rate greater than 0.7MQ/year from BlueTides. Dashed black: star formation rate density of equivalent halos from Illustris 
(|Genel et al.|2014| [Vogelsberg er et al.|2014| ). Solid black: star formation rate density of equivalent halos from MassiveBlack II|Khandai| 
|et al.| ( |201^ Dotted blue: star formation rate density of all halos from BlueTides. Dotted black: star formation rate density of all halos 
from MassvieBlack II. Squares: estimates from HUDF (|Oesch et al.|20I4b|, and HFF A2744( [Oesch et al.|20I4^ . Diamonds: observational 
estimates from CLASH ( Bouwens et al.||20I4b[ [^eng et al.||20I2[ |Coe et al.||20I3| >. Illustris, CLASH, and HUDF lines are reproduced 
from |Oesch et al.| ( |20I4a I. 


while that for MassiveBlack II is a factor of a few larger. 
This suggests that the sub-grid feedback model dominates 
in controlling the star formation rate over both the choice 
of hydrodynamic method and the effect of reionization. It is 
however promising to see, given the differences in the simu¬ 
lations, that the differences are within at most a factor of a 
few in SFR density. 

In Figure we show the cumulative star-format ion rate 
density from galaxies of different halo mass. For comparison, 
we also calculated the cumulative SFR density of Illustris 
at three corresponding redshifts {z = 8,9,10). We see that 
the mass threshold for 50% of star-format ion increases from 
- 5 X at z = 13 to - 4 X at = 8. 

Thus the contribution of small halos to the ionizing photon 
budget becomes increasingly important at higher redshift. 
The larger volume in BlueTides means that it includes ha¬ 
los 10 times more massive than the most massive halos in 
Illustris. These objects contribute less than 10% of the total 
star-format ion density at z = 8,9, and 10. 


4 STELLAR AND AGN UV LUMINOSITY 
FUNCTIONS 

In this Section, we report the stellar AGN luminosity func¬ 
tions in BlueTides. We first describe our source detection 
method, using Source Extractor to validate the results of an 
FoF halo finder at these high redshifts (Section |4.1| ). We then 
describe the stellar UV luminosity function (Section |4. 2.1 ), 
the dust attenuation model (Section [4.2.2| , and the faint-end 



logAAalo [MJh] 

Figure 4. Cumulative star-formation rate density in halos. 
Coloured solid: cumulative SFRD in BlueTides. Cray dashed: 
cumulative SFRD in Illustris at 2 ; = 8, 9,10, from the Illustris 
public data release ( [Nelson et al.||2015| >. Thick black: contour of 
50% cumulative SFRD. Thin black: contour of 10% cumulative 
SFRD. 


slope of the luminosity function (Section [4.2.3| ). We assemble 
the stellar luminosity function from BlueTides and compare 
to observations and other simulations (Section [Olll. F i- 
nally we report the AGN luminosity function (Sect ion |4.3|) . 
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4.1 The identification of Galaxies in BlueTides: 
Source Extractor vs FoF 


The Friend-of-Friends (FoF) algorithm considers only the 
spatial positions of particles and can sometimes artificially 
group dynamically distinct objects into one halo. In this sec¬ 
tion we compare the luminosity function estimated from FoF 
catalogues to that which would be estimated by performing 
standard observational techniques and show that the differ¬ 
ence is small. 

Simulations and observations have long been defining 
objects in different ways. Even the sub-halos in simulations 
do not directly translate to any imaging survey catalogs. In 
imaging surveys, objects are identified by selecting peaks in 
a 2-dimensional image, and the total luminosity depends on 
an aperture radius, (we refer the readers to Stevens et al. 
2014 and references herein) The canonical implementation 
of such an algorithm is Source Extractor (Bertin & Arnouts 


1996); We use SEP, a reimplementation of source extractor 
into python ( Barbary contributors|2014 ). 

We produce a mock 2 — d survey image by projecting 
the BlueTides simulation box along one axis. The star for¬ 
mation rates of particles are distributed into pixels using a 
Gaussian kernel scaled by their SPH smoothing lengths with 


GAEPSI (Peng et al. 2011). We do not attempt to model 


instrumental noise in the mock image. The final star for¬ 
mation surface density image has (2 x 10^)^ pixels, with a 
spatial resolution of 2 h~^Kpc per pixel (~ O.OGarcsec). This 
image is then divided into 100 non-overlapping equal sized 
sub-volumes x — y plane. Each sub-volume has a volume of 
40 X 40 X 400( h“^Mpc)^. This allows us to estimate the cos¬ 
mic variance in e.g. the luminosity functions in fields compa¬ 
rable to those observed. We note that 400/i“^Mpc roughly 
corresponds to a redshift width of Az ~ 1 at z = 8.0, and 
hence each image chunk roughly corresponds to the full vol¬ 
ume of the BoRG survey. 

We run SEP on each of the image chunks, afterwards 
combining them to assemble the full catalogue. Edge effects 
can be safely neglected because the area of the images are 
much larger than the area of the edges. We use a near-zero 
threshold (1 x 10“® /i“^Kpc^) to include all pixels 

that have non-zero star formation surface density. The inte¬ 
grated star formation rate of the objects is measured with 
an aperture size of 8.7 pixels (0.5 arcsec), which roughly cor¬ 
responds to the aperture used by the BoRG survey (Trenti 
et al.|20TT| . Eigure|^ shows visually the process of identify¬ 
ing galaxies using source extractor. We compute UV mag¬ 
nitudes from the star-formation rate of the SEP catalogue 
and the EoE catalogue using ( Stringer et al.||2011 ): 


Muv = -2.5 logio(^) - 18.45 . (2) 

We exclude faint objects with Muv > —14, corresponding to 
unresolved halos typically containing less than 50 dark mat¬ 
ter particles, and for the moment neglect dust extinction, 
which we will discuss in Section [4.2.21 

As shown in Eigure we find that the luminosity 
functions constructed from the EoE halo catalogue and the 
SEP imaging catalogue differ by less than 20%. The main 
differences are that the SEP luminosity function includes 
fewer bright objects than the EoE luminosity function and 
more faint objects. Quantitatively, SEP is ~ 80% of EoE at 





Figure 6. UV luminosity function of galaxies found with source 
extractor compared to that of FoF halos at 2 ; = 8. Green: SEP 
catalogue. Blue: FoF catalogue. Solid grey: the ratio between SEP 
and FoF (axis on right). 


Muv > —22 and ~ 120% of FoF at Muv < —18. These dif¬ 
ferences can be understood by noting that bright (massive) 
FoF halos are often a cluster of smaller objects which SEP 
tends to identify as separate galaxies. The fixed aperture in 
SEP however tends to enhance the UV of smaller objects. 
These effects are smaller at these high redshifts, where large 
groups have not yet formed. Because these differences are 
small, we will use the EoE catalogue for the rest of this 
work. 


4.2 Stellar Luminosity Functions 

4 . 2.1 Intrinsic Stellar Luminosity Functions 

Figure shows the intrinsic galaxy UV luminosity func- 
tionevolving from z = 13 to z = 8. The luminosity function 
is computed from all galaxies within the simulation. Shaded 
areas represent 1 — a uncertainty of the mass function es¬ 
timated from 100 sub-volumes. The luminosity function’s 
evolution in redshift is largely described by an increase in 
amplitude, with minor evolution of the faint-end slope. Fig¬ 
ure also shows a fit of the intrinsic UV luminosity function 
to a modified Schechter model, which captures the simu¬ 
lated luminosity function well at all redshifts. The modified 
Schechter model is 


ln0(M) = In Tin A 

+ A(M"-M)(l-au) (3) 

_ ^q0.1(M*-M) 

This model differs from the usual Schechter model (see 
equation below) only in that the coefficient that deter¬ 
mines the rate of extinction of bright end galaxies is changed 
from 0.4 to 0.1. The best fit values are shown in Tabled The 
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Figure 5. Detecting objects from the mock star-formation intensity image. Background: the star-formation intensity image of the full 
projection of BlueTides (400/i“^Mpc per side). We note that the image is strikingly uniform because of the thickness of the projection. 
Top left: the star-formation intensity image of a single chunk, 40h“^Mpc per side. Top right: all objects identified in a field of view with 
a 10 times zoom, 4h“^Mpc per side. SEP objects are marked in red. Bottom right: a further zoomed-in view of the top right panel, to 
show the identified objects more clearly. 


modified model agrees with the luminosity function at the 
5% level for the whole luminosity range of the simulation. 
Note that the change in the bright end coefficient means that 
one should not directly compare these fit parameters with 
those obtained from a Schechter model. We also do not use 
this model to extrapolate the luminosity function (although 
with a minor 10% diherence in the photon budget, using it 
would not change our conclusions). 


4-2.2 Dust Extinction 


There is evidence that the highest luminosity early galaxies 
are significantly dust obscured ( Wilkins et al.||2Q13 Cen & 
Kimm|2014 ). To produce luminosity functions more compa¬ 
rable to observations, we adopt the screening model from 


Joung, Cen & Bryan (2009). This attenuates the UV lu¬ 


minosity from a pixel in the face-on image of a galaxy by 
a fraction, /uv, proportional to the metal-mass density in 
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SFR [M. ,'yeai] SFR [M.J>-eai] 




Figure 7. Evolution of the intrinsic UV Luminosity function with 
redshift from 2 ; = 8 —13 (colours online). Shaded regions show the 
1 — a sample variance of the mass functions. Dashed lines: best 
fit modified Schechter model. The vertical dash line corresponds 
to the current observation limit of Muv ~ —18. 


Figure 8. Evolution of the UV Luminosity function with dust 
extinction with redshift 2 ; = 8,9,10,11,12,13. (colours online) 
Shaded regions show the 1 — cr uncertainty of the mass functions. 
Dashed lines: best fit modified Schechter model (see Equation]^ 
The vertical dashed line corresponds to the current observational 
detection limit of Mpv ~ —18. 


Table 2. Best-fit parameters of the modified Schechter model for 

the galaxy UV luminosity functions at 2 ; = 8 - 13. Parameters Table 3. Best-fit Schechter Model parameters at 2 ; = 8 - 13 for 

are fit using Equation galaxy stellar UV luminosity functions including dust extinction. 


Parameters are as described in Equation]^ 


z 

(XL 

log 0* 









log 0* 

^uv 

8 

-1.54 ±0.01 

-4.04 ±0.07 

-15.99 ±0.09 

z 

(XL 

9 

-1.59 ±0.02 

-4.17 ±0.14 

-15.44 ±0.19 

8 

-1.84 ±0.03 

-8.90 ±0.21 

-20.95 ±0.15 

10 

-1.55 ±0.04 

-3.66 ±0.18 

-14.09 ±0.28 

9 

-1.94 ±0.03 

-9.74 ±0.24 

-20.77 ±0.17 

11 

-1.51 ±0.07 

-3.42 ±0.21 

-13.03 ±0.40 

10 

-2.01 ±0.04 

-10.27 ±0.33 

-20.39 ±0.22 

12 

-1.40 ±0.07 

-3.42 ±0.10 

-11.78 ±0.37 

11 

-2.07 ±0.05 

-10.77 ±0.42 

-20.00 ±0.27 

13 

-1.32 ±0.10 

-3.82 ±0.08 

-10.98 ±0.44 

12 

-2.12 ±0.06 

-11.40 ±0.50 

-19.65 ±0.31 


13 -2.13 ±0.06 -11.71 ±0.49 -19.11 ±0.30 


that pixel. The value of /uv determines the extinction coeffi¬ 
cient Muv • We apply this dust model to the bright individual 
galaxies at 2: = 8, and find that the dust extinction in UV 
band is fitted by 


Muv - ^uv = exp 


M^v + 22.61 
L72 


(4) 


where Muv is the intrinsic UV luminosity and M^v is the 
UV luminosity with a dust correction. Equation produces 
dust extinction of Muv ~ 1 for Muv = — 21 galaxies, which 
agrees with the upper limit inferred from the UV slope of 


high redshift galaxies by Wilkins et al. (2013) at 2; = 8.0. 


Figure shows the galaxy UV luminosity function with 
dust extinction at 1500A, from 2; = 13 down to 2: = 8. We 
also show the results of a Schechter fit ( Schecht^|1976 ) to 
the observed luminosity functions. The Schechter model is 
widely used to parametrise luminosity functions. We use the 
form provided by |Jaacks et al. (2012). 

ln(/)(M) = ln0* ± InM 

±M(M*-M)(l-au) (5) 


- 10 


0.4(M*-M) 


where M = 0.4 In 10. Parameters are estimated using fit¬ 
ting over In 0, assuming uncorrelated errors (estimated from 


the sub-volumes). The best fit parameters are reported in 
Table [S] Note that the Schechter model does not describe 
the BlueTides luminosity functions at high redshift {z > 10), 
and systematically under-fits the bright end luminosity func¬ 
tion, even after including dust extinction. 


4 . 2.3 Faint-end Slope 

We show the redshift evolution of the faint-end slope of the 
dust extincted galaxy UV luminosity function in BlueTides 
in Figure The best fit model for the evolution is 

aGaia^y{z) = -0.756(1 ± (6) 


The slope of the faint end of the UV luminosity is con¬ 
sistent with that inferred by Bouwens et al. (2014a) and its 
evolution with redshift implies moderate steepening, again 
consistent with an extrapolation of the observed slope evolu¬ 
tion up to 2: = 12, which includes an evolution of M/L ratio 
oc (1±2:)“^‘^ due to the evolution of the Halo Mass Func- 
tior0 The faint end (Muv > —20) of the UV luminosity 


^ in Fig we ignored the 2 : = 10 estimate from [b 


ouwens et 




( 2014a| ), as the authors manually set the faint-end slope at 2 ; = 10 
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many stars in small objects. This is likely due to their stel- 
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Figure 9. Evolution of the faint-end slope of the galaxy UV 
luminosity function. The purple diamonds are the observed slope 
from Hubble legacy surveys ( [Bouwens et al.||2014a| >. Blue circles 
show the slope from BlueTides. 


function in BlueTides is barely affected by the dust extinc¬ 
tion model, thus the redshift-slope relation we give here is 
suitable for inputs of reionization calculations. 


4 . 2.4 Comparison with Observations and Other Models 


Bouwens et al. (2014a) assembled and reanalyzed the UV 


luminosity function evolution from z = 10toz = 4 based 
on all currently available legacy Hubble surveys. The total 
cumulative area is close to 1000 arcmin^, spread over a wide 


redshift range (see also Ellis et al. 2013 Finkelstein et al. 


2014 Trent! et al. 2010). The most recent measurements 
at z = 10 were published in Oesch et al. (2014b). Figure 
compares a compilation of these observational data to 
stellar UV luminosity functions from BlueTides. We show 
the BlueTides luminosity functions extracted from subfields 


with size roughly the area of the BoRG survey (Bouwens 
et al.|20f4a ), the legacy field with the currently largest area. 
As this is a smaller volume than the full simulation, we can 
use the differences between sub-volumes to estimate sample 
variance from current observations, which is shown by the 
shaded areas in FigureThe intrinsic luminosity produces 
more bright galaxies than are observed, a discrepancy which 
is marginally significant compared to cosmic variance. After 
applying a dust extinction correction, this slight tension with 
observations largely disappears, suggesting that dust correc¬ 
tions are indeed significant for the brightest galaxies, even 
at these high redshifts. As we shall discuss in Section |4.3| 
another possibility is that the brightest sources host a sig¬ 
nificant AGN which may make their detection in the galaxy 
samples harder. 

We show comparisons with two other recent high red- 
shift simulations with a maximal volume of 100/i“^Mpc. 


Jaacks et al. ( 2012| ) performed several simulations at differ¬ 
ent resolutions to investigate the shape and slope of high red- 
shift galaxy UV luminosity function. The luminosity func¬ 


tions of Jaacks et al. (2012) have faint end slopes which are 


lar feedback (Ghoi & Nagamine 2011) being insufficiently 


effective at suppressing star formation. We also show the 
corresponding UV luminosity function at 2: = 8,9,10 from 
the public data of Illustris ( Nelson et al.||2015 ). Due to the 
larger volume in BlueTides, Illustris produces fewer bright 
objects and cannot be compared to BlueTides at the most 
massive end. However, at the faint end of the luminosity 
function (Muv < —20), BlueTides and Illustris agree at the 
10% level. This suggests that a stellar feedback model which 
very efficiently suppresses star formation in small halos is 
the most important ingredient when matching the luminos¬ 
ity function at high redshift, just as at low redshift. 


4.3 AGN Luminosity Function 

As described in Section the BlueTides simulation models 
AGN via a self-regulated super-massive black hole model 
following Di Matteo et ah] ( |2005 ). Given a mass accretion 
rate the bolometric luminosity of AGN is 


L = 77 


g^Mbhc^ 

dt 


( 7 ) 


where 77 = 0.1 is the mass-to-light conversion efficiency in 
an accretion disk. 

We convert the bolometric luminosity of AGN in the 


simulation to a UV magnitude using (Fontanot, Gristiani V 
Vanzella|2012 ): 

Tbol 


Muv = —2.51ogiQ 


IbJ^b 


+ 34.1 + Ab,i 


(8) 


where Lbol is the bolometric luminosity of an AGN, fs = 
10.2 ( Elvis et ^|1994 ), and Ab ,uv = —0.48. 

Figure |11 1 shows the UV luminosity function of AGN in 
BlueTides. The luminosity function is cut at Muv = —18.6, 
a limit dictated by the imposed seed mass of our black holes 
(Mseed = 5 X lO^M©). The black hole luminosity function 
is only meaningful for objects that have at least doubled 
their mass since the black hole was seeded, thereby erasing 
the artificially imposed seed mass. For smaller black holes 
the AGN luminosity is significantly suppressed due to the 
artificial absence of black holes in our numerical scheme. 
The AGN luminosity function rises steadily at later times, 
mirroring the evolution in the stellar luminosity function. 
By z > 13 our box contains a negligible number of AGN 
and it is thus impossible to reliably estimate the luminosity 
function. 

Unlike the stellar luminosity function, however, the 
shape of the AGN luminosity function is well-described by 
a power law. We thus fit the AGN luminosity function with 
a power-law model as: 


ln0(M) = ln0* + In A 

+ A(M*-M)(l-au), 


( 9 ) 


where we set a reference magnitude at M* = —18 without 
loss of generality. ( Boyle, Shanks Peterson|1988 Hopkins 


steep compared to observations, producing substantially too 


Richards Hernquist|2007 ) fit the AGN luminosity function 
at low redshifts with a double power law, which allows a bet¬ 
ter fit to the steeper bright end slope. This is not necessary 
for us as the objects which require such a steeper slope are 
brighter than any AGN in BlueTides at 2; = 8.0. The best 
fit parameters of the power law fit are reported in Table 
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log SFR Mg/yr 



Figure 10. Comparing UV Luminosity functions with observations. Solid blue: Intrinsic UV luminosity function in BlueTides. The 
coloured bands shows the cosmic variance in a survey volume of the size of the BoRG survey ( [Bouwens et al.||2014a|>. Sol id red: Dust- 
reddened UV luminosity function in BlueTides. Dotted black: Luminosity functions from the simulations of |Jaacks et al.| ( |201^ . Note 
that they do not provide a simulated luminosity function at 2 ; = 10; also the largest volume at 2 ; = 9 was a 34 /i“^Mpc box. Green wedges: 
intrinsic UV luminosity function from the Illustris public data release [Nelson et al.| ( [20T^ . Black dashed line: intrinsic UV luminosity 
function from the MassiveBlack-II public data jKhandai et al.| ( [2QT4| . Black squares: Observed luminosity functions at 2 : = 8 and 2 : = 10 
from jBouwens et al.j ( [^14a| >. Black diamonds: Observed luminosity function at 2 ; = 9 fromjMcLure et al.j ( [?013| >; [McLeod et al.| ( |2014| . 
Black circles: Observed luminosity function from four bright galaxies at 2 : ~ 10 from [^sch et al.| ( [2014bp 


6.0 6.5 7.0 7.5 B.0 B.5 



Figure 11. AGN UV luminosity function from BlueTides. The 
coloured bands show the 1 — cr sample variance estimated from 
100 sub-volumes (see section[^. We cut the luminosity function at 
the faint end at Muv = —18.6, corresponding to the seed mass of 
the black holes in BlueTides. Dashed lines: the best fit power-law 
model. Symbols: measurements at 2 ; = 5.75 by |Giallongo et al.| 
( [2015[ >; [Kashikawa et al.| ( [2015| >. 


Table 4. Best-fit power-Law parameters at 2 ; = 8 — 12 of the 
AGN UV luminosity function. The fitting model is described by 
Equation [^ 


z 

aL 

log 0* 

8 

-2.45 ±0.02 

-9.33 ±0.04 

9 

-2.36 ±0.07 

-10.44 ±0.09 

10 

-2.35 ±0.04 

-11.91 ±0.07 

11 

-2.02 ±0.24 

-13.79 ±0.33 

12 

-2.43 ±0.22 

-15.30 ±0.28 


5 IMPLICATIONS FOR REIONIZATION 


There are currently few constraints on the process of hydro¬ 
gen reionization. Measurements of the total optical depth 
to the cosmic microwave background (CMB) suggest that 
the redshift of half reionization is 2:haif ~ 10 (Hinshaw et al. 


2013). Small scale CMB experiments have also constrained 


the duration of reionization to be Az < 4.4. IZahn et al.l 
(2012). However, the sources of the UV photons which re¬ 
ionized the universe are subject to extensive debate, with the 
two main candidates being faint galaxies and AGN. Con¬ 
straints on the contribution from different sources can be 
calculated using the measured luminosity functions (See, e.g. 


Fontanot et aL[[2012 

Faucher-Giguere et aL[|2008| [Meiksin 

2005| Haardt & Mac 

au[[2012| Pawlik, Schaye & van Scher- 

penzeel||2009| Bolton & Haehnelt[[2007| [Bunker et aL[[2010| 

Robertson et al.|2013 
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Figure 12. Photon budgets to reionise the universe at 2 ; = 8, as 
a function of the ratio between clumping factor C and escaping 
fraction /esc • Shaded region: UV emissivity from the observed 
luminosity function, extrapolated to Mpv = —10. Black dashed: 
UV emissivity required for a complete reionization at 2 ; = 8. Blue 
solid: UV emissivity from BlueTides, extrapolated to Muv = 
- 10 . 


Quasars have yet to be observed at these high redshifts, 
making the expected impact of AGN on reionization uncer¬ 
tain (see, e.g. Fan et al.|2Q06 Meiksin Madau|1993 Mitra, 


Choudhury Ferrara|20T2 ). However, the most recent con¬ 

straints on the quasar luminosity function from CANDELS 
Goods fields at z ~ 4 — 6 suggest that AGN may make a sig¬ 
nificant contribution to reionization ( Giallongo et al.|2015 ). 
In general, reionization driven by rare bright sources such 
as quasars and large galaxies would progress rapidly, while 
one driven by faint sources would progress more slowly. 

The photon budget can be modelled from first principles 
using the simulated luminosity functions of ionizing sources. 
As ionizing photons may themselves affect the formation of 
small halos ( Madau Pozzetti|2000 ), direct predictions re¬ 
quire simulations which couple radiative transfer and hydro¬ 
dynamics. Furthermore, predicting the UV photon escape 
fraction from galaxies requires exquisite resolution (see, e.g. 
Trac & Gnedin 2011). A simpler approach, which we take 


here, is to use the simulated luminosity functions from a cos¬ 
mological hydrodynamic simulation to estimate the photon 
sources, leaving the escape fraction as a free parameter. 

Our modelling of the photo-ionization rate is based on 


Fontanot et al. (2012), and we refer the reader to this pa¬ 


per for further details of the model (see also Haardt & 
Madau 2012). Here we briefly review the relevant pieces. 


The radiation density Pu[z) of a source species with a evolv¬ 
ing luminosity function (/(Muv,^) and specific luminosity 
Liy(Muv, T-') is 


Pu[z) — / (j)[M\]Y^ z)Lu[M\jy. 

Mcut 


( 10 ) 


For galaxies, the ionizing photo production, Fqal, is based 
on the star formation rate 

Puv^(«)M©year"^Mpc“® 


rGAL(«) = K/es 


1.05 X 1021W • Mpc- 


( 11 ) 


where k = 10 ®^'^ year is the mean opacity, /esc is 

the UV escape fraction. High resolution simulations coupling 
radiative transfer to hydrodynamics in individual galaxies 
have suggested a wide range of possible values for /esc- For 
example, Gnedin, Kravtsov & Ghen (2008) reported that 


high redshift galaxies are very inefficient in releasing their 
photons thus /esc < 5%; Kimm & Gen (2014) suggests 
/esc ~ 0.14. On the other hand, Gen (2005) favours high 


escape fractions, while Yajima, Ghoi & Nagamine (2011) 
suggests that small halos with halo mass < 10^ have 

/esc ~ 0.4. In order to bracket possibilities for the reioniza¬ 
tion contribution from our galaxies we consider two extreme 
scenarios: a low escape fraction model with /esc = 0.05, and 
a high escape fraction model with /esc = 1-0. 

The AGN contribution, Fagn(^), to the ionizing pho¬ 
ton budget depends on the AGN luminosity function and 
associated SED: 


TaGnCz) = / ^—r—^dv, 

hpi' 


( 12 ) 


where = 3.2 x lO^^Hz and z^He = 12.8 x lO^^Hz are the 
ionising frequency of Hydorgen and Helium, hp is the Planck 
constant, and p^^^ is the radiation density of AGN photons. 


Ljy = Ljjy 


(V 

V^UV / 


(13) 


where z/pv = 2 x 10^^Hz is the frequency at 1500A. We will 
adopt ajjv = —1.76 (see, e.g. Hopkins et al.||2007 ). 

The formulae above provide the photon budget. The 
photo-ionization rate required to reionize the universe can 
be estimated theoretically from the recombination rate and 
the clumping factor C{z) as 


Freion( 2() = 0.027z^ — 


c n+z 


7 


V 0.465 J 


(14) 


The clumping factor C{z) describes density variations be¬ 
low the resolution of hydrodynamic simulations and can be 
estimated using higher resolution simulations of the inter- 
galactic medium at the reionization epoch (Einlator et al. 


2012 McQuinn, Oh & Eaucher-Giguere 2011). We use the 
smallest evolving clumping factor from Pawlik et al.| ( |2009| ) 


C{z) = 1 + 432: 


(15) 


Note that a higher clumping factor requires more photons 
to reionize the universe. 

We extrapolate the luminosity functions measured from 
our simulations (see EigurepT| to include the contributions 
from galaxies and AGN smaller than the resolution limit of 
the simulation. We integrate the galaxy luminosity function 
for all UV magnitudes brighter than Muv = —10, the lower 


limit of galaxies that generate ionizing photons (Kuhlen & 
Eaucher-Giguere|2012 ). We consider extrapolating the AGN 
UV luminosity to Muv = —12 from Muv = —IS, the 
faintest AGN in the BlueTides simulation. However, further 
decreasing the AGN threshold does not increase the number 
of ionizing photons significantly. 

In Eigure we show the luminosity density in Blue¬ 
Tides at 2: = 8 and the photon budget to reionize the uni¬ 
verse at 2: = 8 as a function of the ratio between the 
clumping factor and escaping fraction. The UV luminosity 












































































12 Yu Feng et. al 


density in BlueTides, like the luminosity function, is con¬ 
sistent with observations. We also see that if galaxies alone 
reionise the universe by z = 8, BlueTides implies a ratio of 
= 4. Assuming a clumping factor of Ciz = 8) = 2, 
reionization at z = 8 requires a high escape fraction of 
/esc ~ 60%. With less UV photon production, it would be 
difficult for the galaxies in BlueTides alone to reionize by 


Recent observations (especially Planck) favour reioniza¬ 
tion at z < 8. Bouwens et al. (2015) analyzed the constraints 
on the ionization history, incorporating recent Planck results 
with constraints from quasar absorption and Lyman-a emis¬ 
sion line measurements. Figure shows a comparison be¬ 
tween the ionization rate in BlueTides and the observational 
constraints under various scenarios for the AGN luminosity 
and galaxy escape fraction. 

As mentioned above, provided /esc = 0.5, galaxies alone 
can produce an ionization history which completes by z = 
8.0, consistent with current observational constraints. With 
a smaller UV escape fraction, the contribution from faint 
AGN instead drives most of reionization, dominating over 
the UV photons from galaxies. As seen in the upper panels 
of Figure AGN tend to produce an ionization rate which 
increases faster than that from galaxies. It is also interesting 
to note that the AGN scenario implies that faint AGN exist 
with luminosities down to Muv = —12, corresponding to a 
blackhole mass of 10^ M© (assuming Eddington accretion), 
and much smaller than is currently observed. Overall, it is 
difficult to produce reionization completing by z = 8 without 
either invoking a high UV escape fraction from galaxies or 
a significant contribution from a faint AGN population. 


6 CONCLUSIONS 

We have performed BlueTides, a high resolution, 
400^/i“^Mpc^ uniform volume hydrodynamical simu¬ 
lation. BlueTides includes a pressure-entropy formulation 
of smoothed particle hydrodynamics, gas cooling, star 
formation (including molecular hydrogen), black hole 
growth and models for stellar and AGN feedback processes. 
BlueTides is the first cosmological large volume hydro 
simulation to incorporate a “patchy “ reionization model 
producing an extended hydrogen reionization history. We 
have reported the high redshift (z > 8) UV luminosity 
functions of galaxies and AGN in BlueTides, and examined 
the implications for reionization. 

We find good agreement between the expected star for¬ 
mation rate density in BlueTides and current observations 
at 8 < z < 10. By using the star formation rate in our galax¬ 
ies we make predictions for the intrinsic galaxy luminosity 
functions and show that they compare favourably to obser¬ 
vations from Hubble Space Telescope (HST) legacy fields. 
The brightest galaxies are yet to be observed and we pre¬ 
dict that upcoming larger area surveys should start detect¬ 
ing them. At z = 8 some dust may be required to reproduce 
the currently observed bright end in the HST surveys. 

Our simulation predicts a faint-end slope of the lumi¬ 
nosity function consistent with observations. When fit to 
a Schechter luminosity function, the slope varies between 
1 .8atz = 8 tO(a~ —2.1 at z = 10 with an evolution 
in the slope oc (1 + . The AGN luminosity functions 


from BlueTides can be fit by a power law with a slope con¬ 
sistent with the most recent observations from GANDELS 
Goods fields (Giallongo et al. 2015). The AGN population 
evolves quickly at these redshifts with the brightest quasars 
reaching Muv ~ — 25 at z ~ 8, which is at least an order of 
magnitude fainter than SDSS quasars at z ~ 6. By combin¬ 
ing the AGN and galaxy luminosity functions we find that 
the bright end Muv ~ —21.5 flattens due to the increased 
AGN activity above Muv ~ —21.5. 

We find that a high (~ 50%) escape fraction is still 
required for galaxies alone to produce enough photons to 
reionize the Universe by z = 8. Our high escape frac¬ 
tion model supports the conditions proposed by |Kuhlen 


Eaucher-Giguere (2012): a reionization model that includes 


mostly galaxies requires both an extrapolation to very faint- 
end Muv = —10 and a sharp increase of the escape fraction 
(up to 50%) at high redshift, in agreement with [Bouwens 


jet al.j ( |2015| ). Eor lower escape fractions (closer to 10%- 
20%), a possible source of the extra photons are fai nt AGN 
with black hole masses as faint as M ~ lO^M© (Madau 
jet al.|[20Q4 ). Giallongo et al. (2015) has suggested that the 
faint end of the AGN luminosity function at z < 5.75 may 
favour these sources. AGN would lead to a relatively quick 
reionization which could soon be testable observationally. 
Alternatively, reionization may not complete until < 8, as 
suggested by Einkelstein et al. (2014). and consistent with 


the recent optical depth measurements from Planck (Planck 
Gollaboration et al.||2015 ). 
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